The  Laplace  Transform  is  frequently  encountered  in  mathematics,  physics,  engineering  and 
other  fields.  However,  the  spectral  properties  of  the  Laplace  Transform  tend  to  complicate 
its  numerical  treatment;  therefore,  the  closely  related  “Truncated”  Laplace  Transforms  are 
often  used  in  applications. 

We  have  constructed  efficient  algorithms  for  the  evaluation  of  the  Singular  Value  Decomposi¬ 
tion  (SVD)  of  Truncated  Laplace  Transforms;  in  the  current  paper,  we  introduce  algorithms 
for  the  evaluation  of  the  right  singular  functions  and  singular  values  of  Truncated  Laplace 
Transforms.  Algorithms  for  the  computation  of  the  left  singular  functions  will  be  introduced 
separately  in  an  upcoming  paper. 

The  resulting  algorithms  are  applicable  to  all  environments  likely  to  be  encountered  in 
applications,  including  the  evaluation  of  singular  functions  corresponding  to  extremely  small 
singular  values  (e.g.  10~1000). 
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1  Introduction 


The  Laplace  Transform  £  is  a  linear  mapping  L2[0,  oo)  — >  L2[ 0,  oo);  for  a  function  /  G  L2[0,  oo), 
it  is  defined  by  the  formula: 

COO 

(£(/))  M  =  /  e-to/(t) df.  (1) 

Jo 

As  is  well-known,  £  has  a  continuous  spectrum,  and  £_1  is  not  continuous  (see,  for  example, 
[lj).  These  and  related  properties  tend  to  complicate  the  numerical  treatment  of  £. 

In  addressing  these  problems,  we  find  it  useful  to  draw  an  analogy  between  the  numerical 
treatment  of  the  Laplace  Transform,  and  the  numerical  treatment  of  the  Fourier  Transform  T\ 
for  a  function  /  G  L1(M),  the  latter  is  defined  by  the  formula: 

/oo 

(2) 

-OO 


where  w  G  M. 

In  various  applications  in  mathematics  and  engineering,  it  is  useful  to  define  the  “Trun¬ 
cated”  Fourier  Transform  Tc  :  L2[— 1,1]  — >•  £2[— 1,1];  for  a  given  c  >  0,  Tc  of  a  function 
/  G  £2[— 1,  1]  is  defined  by  the  formula: 

ow))  m = yV""/(<)d(.  (3) 

The  operator  Tc  has  been  analyzed  extensively;  one  of  most  notable  observations,  made  by 
Slepian  et  al.  around  1960,  was  that  the  integral  operator  Tc  commutes  with  a  second  or¬ 
der  differential  operator  (see  |2j).  This  property  of  Tc  was  used  in  analytical  and  numerical 
investigations  of  the  eigendecomposition  of  this  operator,  for  example  in  [21 0 1  El  El  U\  IS]. 

For  0  <  a  <  b  <  oo,  the  linear  mapping  Ca ^  ■  L2[a ,  b]  — >  L2[ 0,  oo)  defined  by  the  formula 

(■ £a,b(f ))  M  =  [  e~tuJf(t)dt,  (4) 

J  a 

will  be  referred  to  as  the  Truncated  Laplace  Transform  of  /;  obviously,  Ca,b  is  a  compact 
operator  (see,  for  example,  DP)  - 

The  Singular  Value  Decomposition  (SVD)  of  Ca,b  has  been  analyzed,  inter  alia,  in  [1]  and 
f9|;  Bertero  and  Griinbaum  observed  that  each  of  the  symmetric  operators  (£a,b)*  °  £a,&  &nd 
-Ca,ft°(£a,b)*  commutes  with  a  differential  operator  (see  |9j).  Despite  [QlfTOlUl  11,  12,  1 3  11  1. '15j, 
much  more  is  known  about  the  numerical  and  analytical  properties  of  Tc  than  about  the 
properties  of  £a,&- 

We  have  constructed  algorithms  for  the  efficient  evaluation  of  the  of  the  SVD  of  Ca^,-  In  this 
paper,  we  introduce  algorithms  for  the  efficient  evaluation  of  the  right  singular  functions  and 
singular  values  of  £a,&-  The  remaining  algorithms,  including  the  algorithm  for  the  numerical 


1 


evaluation  of  the  left  singular  functions,  will  be  discussed  in  upcoming  papers  along  with 
additional  analytical  results. 

The  paper  is  organized  as  follows.  Section [2] summarizes  the  various  standard  mathematical 
facts  and  simple  derivations  that  are  used  later  in  the  paper.  Section[3]contains  the  derivation  of 
various  properties  of  the  right  singular  functions  of  the  Truncated  Laplace  Transform,  which  are 
used  in  the  algorithms.  Section[4]describes  the  algorithms  for  the  evaluation  of  the  right  singular 
functions  and  singular  values  of  the  Truncated  Laplace  Transform.  Section [5] contains  numerical 
results  obtained  using  the  algorithms.  Section  [6]  contains  generalizations  and  conclusions. 
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2  Preliminaries 

2.1  The  Legendre  Polynomials 

In  this  subsection  we  summarize  some  of  the  properties  of  the  the  standard  Legendre  Polynomi¬ 
als,  and  restate  these  properties  for  shifted  and  normalized  forms  of  the  Legendre  Polynomials. 

We  define  the  Shifted  Legendre  Polynomial  of  degree  k  =  0, 1, which  we  will  be  denoting 
by  Pf,  by  the  formula 

P*k(x)  =  Pk(  2x-l),  (5) 

where  Pk  is  the  Legendre  Polynomial  of  degree  k\  the  standard  definition  of  the  Legendre 
Polynomials  can  be  found,  inter  alia,  in  fl6l. 

As  is  well-known,  the  Legendre  Polynomials  form  an  orthogonal  basis  in  L2[— 1, 1],  but 
they  are  not  normalized;  it  immediately  follows  that  the  Shifted  Legendre  Polynomials  form 
an  orthogonal  basis  in  L2[0, 1]  and  that  they  are  also  not  normalized.  Therefore,  we  find  it 
convenient  to  define  the  Normalized  Shifted  Legendre  Polynomial  of  degree  k  =  0, 1, ...,  which 
we  will  be  denoting  by  Pf.  by  the  formula 

P*(x)  =  Pf(x)V2k  +  l-  (6) 

the  Normalized  Shifted  Legendre  Polynomials  Pq  ,  Pf, ...  form  an  orthonormal  basis  in  L2[0, 1]. 

The  following  well-known  properties  of  the  Legendre  Polynomials  can  be  found,  inter  alia, 
in  HS|,  PH: 


( Pfe(x))2dx 


2 

2k  +  1 


(7) 


(k  +  l)Pk+i(x)  =  (2k  +  l)xPk(x)  -  kPk-i(x) 


(8) 


JL  d 


(1  -  x^)—Pk(x)  =  —kxPk(x)  +  kPk_i(x ) 
dx 


(9) 


d  (  ,.  2.  d 


—  (  (1  -  x2)— Pk(x)  )  =  -k(  1  +  k)Pk(x) 


(10) 


(2k  +  1  )Pk(x)  =  —  (Pk+ i(x)  -  Pfc_i(x)) 


(11) 
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The  following  properties  of  the  Shifted  Legendre  Polynomials  are  easily  derived  from  the 
properties  of  the  Legendre  Polynomials  by  substituting  ([5])  into  ([7  11). 


(  KW)Idl  =  iTT 


(12) 


1  fkP*(x)  (1  +  k)P*(x) 

^-2  lW  +  ^)+ 


(13) 


x(1 "  x)TxP*k{x)  =  WTW)  (PU{x)  ~  Pk+l{x^ 


(14) 


^  (x(l  -  x)  Ap*(x)  )  =  _fc(  1  +  fc)Pfc*(s) 


(15) 


2.2  The  Legendre  Functions  of  the  second  kind 

As  is  well-known,  the  Legendre  Polynomial  Pk{x )  is  not  the  only  solution  for  the  differential 


equation  (10)  in  the  interval  [—1,1];  the  other  solution  is  the  Legendre  Function  of  the  second 
kind  Qk(x),  defined  by  the  formula 


Qk(x)  =  7}  J  (x-t)  1Pk(t)dt, 


(16) 


where  Pk  is  the  Legendre  Polynomial. 

Having  defined  the  Shifted  Legendre  Polynomials,  we  find  it  convenient  to  similarly  define 
the  Shifted  Legendre  Function  of  the  second  kind  of  degree  k,  which  we  will  be  denoting  by  Q*k , 
by  the  formula 


Ql(x)  =  Qk{ 2®  -  1). 

The  following  identities  can  be  found,  for  example,  in  m,  m 

Qk(z)  =  (-1  )k+1Qk(-z), 


(17) 


(18) 


d(p 


roo 

Qk(z)  =  /  - - 

Jo  [z  + 

By  Q,  ([16]),  ([18])  and  (]T7|), 


COS 


k-\- 1  ' 


(19) 


(x  +  y)~1P*{x)dx  =  2(-l)kQl{y  +  l)V2k  +  l 


(20) 


for  all  y  >  0. 
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2.3  Singular  Value  Decomposition  (SVD)  of  integral  operators 

The  Singular  Value  Decomposition  (SVD)  of  integral  operators  and  its  key  properties  are 
summarized  in  the  following  theorem,  which  can  be  found,  for  example,  in  |18|. 

Theorem  2.1.  Suppose  that  the  function  K  :  [c,d\  x  [a,b\  — >  M  is  square  integrable,  and 
T  :  L2[a,  b]  — >  L2[c ,  d]  is  defined  by  the  formula 

(T(/))0)=/  K(x,t)f{t)dt.  (21) 

J  a 

Then,  there  exist  two  orthonormal  sequences  of  functions  uo,ui,...,  where  un  :  [a,  6]  — >  M  and 
vo,vi,...,  where  vn  :  [c,d]  — >  M,  and  a  sequence  ao,ai,...  £  M,  where  chq  >  a\  >  ...  >  0,  such 
that 


(T(/))  (x)  =  y  ,an 


OO 

E' 

71= 0 


un(t)f(t)dt  vn(x 


(22) 


for  any  f  £  L2[a,b\.  The  sequence  otQ,ot\,...  is  uniquely  determined  by  K. 


The  functions  uq,  u±,  ...  are  referred  to  as  the  right  singular  functions,  the  functions  vo,v\, ... 
are  referred  to  as  the  left  singular  functions,  and  the  values  o.q ,a\,...  are  referred  to  as  the 
singular  values  of  the  operator  T.  Together,  the  right  singular  functions,  the  left  singular 
functions  and  the  singular  values  are  referred  to  as  the  SVD  of  the  operator  T. 


It  immediately  follows  from  Theorem  2.1  that 


TiyUn)  —  anvn, 


(23) 


T*(vn)  =  anun. 


(24) 


Observation  2.2.  The  right  singular  functions  uq,  u\,  ..  of  T  are  eigenfunctions  of  the  operator 
T*  o  T  and  the  left  singular  functions  v$,v\,...  are  eigenfunctions  of  the  operator  T  o  T*;  the 
singular  values  ao,  a\, ...  of  T  are  the  square  roots  of  the  eigenvalues  of  T*  o  T  and  T  oT* .  In 
other  words,  for  every  n  =  0, 1, ..., 


and 


((foT)y(r)  =  /  K(x,t)[  /  K{x,t)un{t)dt)  dx  =  a2nun{r)  (25) 


((T  o  T*)(yn))  (£)  =  J  K(£,t)  K{x,t)vn(x)dxSj  dt  =  a2nvn{£)  (26) 
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Remark  2.3.  The  function  K  can  be  expressed  using  the  singular  functions  as  follows  (see 

w>, 


oo 

K(x,t)  =  y ^vn(x)anun(t)  (27) 

n= 0 

and  it  can  be  approximated  by  truncation  of  small  singular  values  (also  see  |18|): 

p 

K(x,t )  ~  y ^vn(x)anun(t)  (28) 

n=0 
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3  Analytical  apparatus 

3.1  Bounds  on  the  Legendre  functions  of  the  second  kind 

For  any  x  >  1,  the  function  Qk{x)  (defined  in  ([l6|))  decays  rapidly  as  k  grows.  More  formally, 
for  any  6  >  0  and  k  =  0, 1, there  is  a  uniform  bound  on  \Qk(x)\,  where  x  >  1  +  <5;  the  bound 
decreases  super  algebraically  as  k  grows.  The  following  lemma  provides  an  explicit  bound. 

Lemma  3.1.  Suppose  that  5  >  0;  then,  for  all  y  >  0, 


|Qfc(l  +  d  +  y)\  <  \  log  (  2 — ~ —  ]  +  1 


1  +  5 


k-\- 1 


where 

5  =  ^(l  +  5)2 -1 
and  Qk  is  defined  in  0- 


Proof.  By  (16), 


\Qkifi  +  &  +  y)\  —  I 
■Jo 


dap 


0  ^(1  +  S  +  y)  +  cosh(y?) y/(l  +  8  +  y)2  -  lj 


Since  (1  +  <5  +  y)  >  (1  +  5), 

r 

|Qfc(i  +  d  +  y)\  < 

■Jo 


d(p 


((1  +  <5)  +  cosh(<^)y/(l  +  5 )2  —  1^ 


fc+i 1 


Since  6  >  0,  clearly  5  >  0  and  therefore, 


\QkO-  +  d  +  y)\  < 


roo 

Jo  ( L 


dp 


+  d  cos. 


k+ 1  ' 


Introducing  the  notation 


v  =  log  2 


l  +  <5 


(29) 


(30) 


\  k+ 1  ‘ 


(31) 


(32) 


(33) 


(34) 


we  break  the  integral  in  (33)  into  integrals  on  the  two  intervals  [0,  &*)  and  [v,  oo): 

dp  dp 


\Qk(fi  +  d  +  y)\  < 


+ 


\  k-\- 1  ' 


0  ^  1  +  6  cosh(</?)^  "  Jv  |^1  +  5cosh((/s)J 


f 


(35) 
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Clearly, 


and 


so  that, 


< 


k-\- 1  —  /  k-\- 1  ! 

1  +  5  coshfo?) )  ( 1  +  5  , 


fc+i  — 

1  +  5  cosh((p) )  (  5  ex 


k+ 1  ' 


IQfc(i  +  5  +  y)l  < 


1  +  5 


k-\- 1 


+ 


f 


d  ip 


5  ex 


k+ 1 ' 


Substituting  (34)  into  (38),  we  obtain 


\Qkif  +  5  +  y)\  < 


1 


1  +  5 


fe+i 


log  ( 2— i—  |  +  ' 


k  +  i  r 


and  from  it,  we  obtain  (29). 


Corollary  3.2.  By,  (Tty  and  (29) 


I Qfe(l  +  5/2  +  y)|  <  (  log  (  2  ]  +  1 


1  +  5 


k-\- 1 


where 

5=  v/(l  +  5)2-l, 
5,  y  >  0  and  Q*k  is  defined  in  m- 


(36) 

(37) 

(38) 

(39) 

□ 

(40) 

(41) 


3.2  The  Truncated  Laplace  Transform 

Definition  3.3.  For  any  pair  of  real  numbers  a,b,  such  that  0  <  a  <  b  <  oo,  the  Truncated 
Laplace  Transform  La g  is  the  linear  mapping  L2[a,b\  — >  L2[0,oo),  defined  by  the  formula 

(£,,&(/))  M  =  fe-fm,  (42) 

J  a 

Obviously,  the  adjoint  of  Lag,  is 

roc 

(( La,bY(g))(t)=  e-^g^du.  (43) 

Jo 

The  operators  Lag  and  (La,b)*  are  compact,  the  range  of  (La,b)*  is  dense  in  L2[a,  b]  and  the 
range  of  Lag  is  dense  in  L2[0,  oo)  (see,  for  example,  [I]). 


3.3  The  SVD  of  the  Truncated  Laplace  Transform 


By  Theorem  2.1  there  exist  an  orthonormal  sequence  of  right  singular  functions  uo,ui,...  E 
L2[a,  b],  an  orthonormal  sequence  of  left  singular  functions  vo,v\, ...  E  L2[ 0,  oo)  and  a  sequence 
of  real  numbers  aq, ...  E  M  such  that 


°°  /  rb  \ 

(£a,b(f))  M  =  X]  (  /  un(t)f(t)dt  I  Vn(uj), 
n= 0  ' 


and  for  all  n  =  0,1, 


(44) 


(45) 


and 


(£ a,b )  (^n)  —  ^n^ni 


(46) 


CX-n  ^  Ckn+1  ^  0.  (47) 

Remark  3.4.  The  multiplicity  of  the  singular  values  of  £a  &  is  one  (see  IS]);  in  other  words, 
for  all  n  =  0,1, ... 


OLn  '5>  ®n+ 1- 


(48) 


Remark  3.5.  According  to  Observation  2.2,  the  right  singular  functions  uo,ui, ...  of  Ca,b  are 
eigenfunctions  of  the  integral  operator  (£a,b)*  °  £-a.b  '■  L2[a ,  b ]  — >  L2[a,  6]  given  by  the  formula 


Wa,bT  °  £«fi)  (/))  (t)  =  l°°  e-*  (J  e-“sf(s)d s^j  d^  =  j"  ^rsmds, 

(49) 

and  the  corresponding  eigenvalues  of  (Ca^)*  o  Ca  b  are  a^,  a.2, ...,  where  an  is  the  singular  value 
of  Ca,b  associated  with  the  right  singular  function  un.  In  other  words, 

(((£a,b)*  °  £a,b)  (un))  (t)  =  [  -^—Un{s)ds  =  a2nlin{t) .  (50) 

Ja  t  +  S 

Similarly,  the  left  singular  functions  vn  of  Ca  b  are  eigenfunctions  of  the  integral  operator 
£a,b  °  ( £a,b )*  '■  L2[ 0,  oo)  — >  L2[ 0,  oo)  given  by  the  formula 


{{Cafio  {Ca,b)*)  (g))  (w)  = 


=  /  e 


—U)t 


1  \  r 

e~ptg(p)dpj  dt  =  j 


oo  —a(u>+p)  _  -b(<jj+p) 


-g(p)dp, 


(51) 


/o  /  Jo  u  +  P 

and  the  corresponding  eigenvalues  Ca ^  o  (£ab)*  are  otQ,  ai>  ••••  hi  other  words, 

fOO  —a(w+p)  _  e—b{w+p) 

/  0  OJ  +  p 


((Ca,b  o  (Ca,b)*)  (vn))  (w)  =  / 

Jo 


vn(p)dp  =  ot2nvn(u).  (52) 
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3.4  The  differential  operators  Dt  and  Du  associated  with  the  singular  func¬ 
tions  of  Ca<b 

In  this  subsection  we  summarize  several  properties  related  to  the  differential  operator  Dt , 
defined  by  the  formula 

(a(/))  (i)  =  ^  (V  -  a2)(b2  -  *2)^/(<)^  -  2(^2  -  a2)f{t),  (53) 

where  /  E  C2[a,b\ ;  and  properties  related  to  the  differential  operator  defined  by  the 
formula 


D„ 


H  = 


_cP_ 

do;2 


_cP_ 

did2 


f(u)J  +  (a2  +  62)-j^  +  (-a2b2w2  +  2a2)  /(a;), 

(54) 

where  /  E  C4[0,  oo)  n  L2[0,oo).  For  a  derivation  of  these  properties,  see  [9]. 

Theorem  3.6.  The  differential  operator  Dt,  defined  in  (53),  commutes  with  the  integral  oper¬ 
ator  (£a,b)*  °  £a,b>  (specified  in  (49))  in  L2[a,b\.  In  other  words, 

Dt  o  ((£a,b)*  °  Ca,b)  =  (( £a,b )*  °  £a,b )  °  A  (55) 


Theorem  3.7.  The  differential  operator  Du,  defined  in  (54),  commutes  with  the  integral  op¬ 
erator  £a,b°  (£a,b)* ;  (specified  in  (51))  in  L2[0,oo).  In  other  words, 


£a,b  °  (£a,b)*  °  4  =  4°  A, b  °  ( £<x,b T 


(56) 


Theorem  3.8.  The  right  singular  functions  uq,u\,  ...  (defined  in  (45))  of  £a_b  (defined  in  (. 4 2)) 
are  also  the  eigenfunctions  of  Dt. 


Theorem  3.9.  The  left  singular  functions  vq,vi,...  (defined  in  (46))  of  £a^  (defined  in  (4%)) 
are  also  the  eigenfunctions  of  Du . 

We  denote  the  eigenvalues  of  the  differential  operator  Df  by  Xo,Xi,  ■■■,  and  the  eigenvalues 
of  the  differential  operator  Du  by  XoiXii  ■■■■  By  Theorem  3.8,  the  singular  function  un  is  the 
solution  to  the  differential  equation 


dt 


tt  (  (i2  -  a2)(b2  -  t2)D-un(t)  )  -  2 (f2  -  a2)un(t )  =  xnun(t), 


dt 


(57) 


and  by  Theorem  3.9,  the  left  singular  function  vn  is  the  solution  to  the  differential  equation 


d2  /  2  d2 


do;2  \  do;2 
=  X*kVk(u). 


dw 


UJ- -—vk(uj)  +  (a2  +  b-)—  a r—Vk(uj)  +  {-a~bzu2  +  2a~)vk(u)  = 


da; 


(58) 
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Remark  3.10.  The  singular  values  an  (defined  in  (45))  of  the  integral  operator  Ca,b  are 
known  to  decay  exponentially  as  n  grows;  consequently,  the  direct  numerical  computation  of 
the  singular  functions  of  Ca,b  beyond  the  first  few  singular  functions  is  impossible. 

The  differential  operators  Dt  and  Du  are  advantageous  in  the  numerical  treatment  of  the 
singular  functions  un  and  vn  because  their  eigenvalues  increase  with  n,  and  because  the  differ¬ 
ential  operators  can  be  treated  using  numerical  tools  developed  for  differential  equations.  Such 


tools  are  developed  below  in  Sections  3A 
of  the  operator  La^  in  Section  |4j 


3.5  3.6  3.71  and  3.8  and  used  to  construct  the  SVD 


3.5  The  operator  T1  and  the  function  77 


The  right  singular  functions  un  (see  (45))  of  Ca (see  (42))  are  defined  of  the  interval  [a,  &];  we 
find  it  convenient  to  shift  this  interval  to  the  interval  [0, 1]. 

We  introduce  the  operator  T7  :  [0, 1]  — >  [0,  oo),  defined  be  the  formula 


(T7(/))(£>)=  f1  e  *(x+i-i)f(x)dx. 
Jo 

This  operator  is  related  to  the  operator  where 
7  =  b/a 


(59) 


(60) 


by  a  the  changes  of  variables 


t  —  a 
b  —  a 


t  =  a  +  (b  —  a)x, 


and 


Cj  =  u)(b  —  a)  , 


uj  = 


u 

b  —  a 


(61) 


(62) 


We  denote  the  singular  values  of  T7  by  do >07,...,  the  right  singular  functions  of  T7  by 

i/jq,  tpi, ...  and  the  left  singular  functions  of  T7  by  vq,vi, . 

Suppose  now  that  0  <  a  <  b  <  00.  Then  a  simple  calculation  shows  that  for  any  n  = 

0,1,2,... 


ipn(x)  =  Vb  —  a  un(a  +  (b  -  a)x), 

(63) 

vn(uj)  =  - vn(u/(b  a)), 

yo  —  a 

(64) 

(up  to  the  ambiguity  in  sign)  and 

(65) 
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where  an,un  and  vn  are  a  singular  value,  right  singular  function  and  left  singular  function  of 
Ca^,  and  where  an ,ipn  and  vn  are  a  singular  value,  right  singular  function  and  left  singular 
function  of  T7  and  7  =  6/0. 

The  operator  T*  o  T7  is  defined  by  the  formula 


(  (T*  O  Try)  (/))  (x)  = 


1 


x  +  y  +  (3 


f(y)dy, 


(66) 


with  /  E  L2\ 0, 1],  and  j3  is  defined  by  the  the  formula: 


P  = 


2  a 


7  —  1  b  —  a 


(67) 


By  Observation  2.2  the  right  singular  functions  ip o,ipi, ...  of  T7  are  the  eigenfunctions  of  the 
integral  operator  T*  o  T7.  Clearly,  T*  o  T7  has  the  same  eigenvalues  as  (Ca,b)*  o  Ca^: 


( (T*  o  T7)  (ipn))  (x)  = 


1 


x  +  y  +  f3 


ipn(y)dy  =  a2nipn(x). 


(68) 


Similarly,  by  (53)  and  (61),  ipo,ipi, ...  are  the  eigenfunctions  of  the  differential  operator  Dx, 
defined  by  the  formula 

(Dx (/))  (x)  =  ^  (x(l  -  x)(/3  +  x)(/3  +  1  +  x)^;/(;r)l  -  2x(x  +  fi)f(x)-  (69) 


In  other  words,  ipn  is  the  solution  to  the  differential  equation 


_d_ 

dx 


x)(/3  +  x)(/3  +  1  +  x) 


2x(x  +  (3)lpn  (x)  -  XnV'n(x)  =  0, 

(70) 


with  XOj  Xi>  the  eigenvalues  of  Dx. 

Remark  3.11.  A  simple  computation  shows  that  the  eigenvalues  xo5  Xi ,  of  the  operator  Dx 
are  related  to  the  eigenvalues  X0)Xi5  •••  of  the  operator  Dt ,  defined  in  (53)  by  the  formula 


Xn  =  {b-a)2Xn- 


(71) 


Remark  3.12.  The  operator  T7  is  determined  by  the  single  parameter  7;  therefore,  the  singular 
value  decomposition  of  the  operator  Ca^a  is  determined  by  7,  in  the  following  sense.  The 
sequence  of  singular  values  07,  ai, ...  of  the  truncated  Laplace  transform  La^a  depends  only  on 
7,  and  is  independent  of  the  value  of  a.  If  7  =  ^  =  |,  then  the  the  sequences  of  right  and  left 
singular  functions  of  Ca ^  are  identical  to  those  of  the  truncated  Laplace  transform  C~  r  up  to 
trivial  scaling. 
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3.6  Expansion  of  fi>n  in  the  basis  of  Legendre  Polynomials 

Let  /  be  a  smooth  function  in  L2[0, 1],  then  f  can  be  expressed  in  the  basis  of  Normalized 
Shifted  Legendre  Polynomials  P£  (defined  in  (p|);  let  h  =  (ho,  hi,  ...)T  be  the  vector  where 

hk=  f  f(x)P£(x) dx,  (72) 

Jo 

then  clearly  h  is  the  vector  of  coefficients  in  the  expansion, 

OO 

/(*)  =^2hkPk(x)-  (73) 

fc=0 


We  introduce  the  notation  hn  =  (/iq  ,  h™,  ...)T  for  the  vector  of  coefficients  of  the  expansion 
of  the  function  (defined  in  (|63|) )  in  the  basis  of  Normalized  Shifted  Legendre  Polynomials; 
where  the  element  h ^  is  dehned  by  the  formula 

K=  [  ipn(x)P£(x)dx,  (74) 

Jo 

so  that 


^n(x)  =  J2hkPk(x)- 

k= 0 


(75) 


3.7  Decay  of  the  coefficients 


Since  the  function  (defined  in  (63))  is  a  smooth  solution  of  a  differential  equation  (specified 
in  @0,  we  expect  the  coefficients  /ij)  in  the  expansion  of  ipn  to  decay  rapidly.  In  this  subsection 
we  provide  an  estimate  for  the  actual  decay. 

Lemma  3.13.  Suppose  that  0  <  /3  <  oo  .  Then, 


1 


o  o  x  +  y  +  fJ 


P£(x)dx  )  dy  < 


2i/2fc  +  1 
fc+i 

1  + 


log  2 


where  P£  is  defined  in  and 

P  =  \](l  +  m)2-l  =  2y/Ml+p). 


l+l 

p  . 


+  1 


(76) 


(77) 
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Proof.  Based  on  (20) 


(x  +  y  +  P)  Pt{x)  dx 


—  2  Q(,(y  +  P  +  l)y/2n  +  1, 


where  Q*k  is  defined  in  ([lTj) .  So,  by  Corollary  3.2 


(■ x  +  y  +  P )  1Pk(x)dx 


< 


2V2k  +  1 

fc+i 

1  + 


log  2 


1  +  P 

P  . 


+  1 


By  squaring  (79)  and  integrating  over  y,  we  obtain  (76). 


(78) 


(79) 


□ 


Lemma  3.14.  Suppose  that  hk  is  the  k  +  1-th  coefficient  in  the  expansion  of  ipn,  specified  in 
(74);  then, 


\hl\  <  cnyj2k  +  1 


1 


1  +  P 


fc+i 


where,  cn  is  defined  by  the  formula 

1  +  P 


cn  =  2an2  log  (  2- 


and  P  is  defined  by  the  formula 


P 


+  1  —  2an  2  log  2  + 


I  y  —  1 

2(7  +  1) 


P  =  V(i  +  rn2)  - 1  =  7)  =  2V2^I. 


(80) 


+  1  ,  (81) 


(82) 


The  parameters  an,  P  and  7  in  (81)  and  (82)  are  defined  in  (45),  (61)  and  (60),  respectively. 


Proof.  We  substitute  (68)  into  (74)  and  change  the  order  of  integration: 

-1  /•! 


tit  =  a  2 


1 


0  Jo  x  +  y  +  P 


Pk(x)'tpn(y)dxdy  = 


=  an2  Jo  My)  „  ,  ].  ,  aPk(x)dx  )  dy. 


x  +  y  +  P 


(83) 


By  the  Cauchy-Schwarz  inequality, 

\K\<an2\  [  (Mv))2 dy\ 


1  /  /•! 


1 


'0  x  +  y  +  P 


Pf(x)dx)  dy. 


(84) 


Now,  by  (76), 


\hk\  <  an2Vl 


( 


\ 


2V2  k  +  1 

fc+i 

1  + 


log  2 


1  +  P 


+  1 


(85) 


□ 
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3.8  A  matrix  representation  of  the  differential  operator  Dx  in  the  basis  of 

>*" 


The  purpose  of  this  subsection  is  to  express  the  differential  operator  Dx  in  the  basis  of  Nor¬ 
malized  Shifted  Legendre  Polynomials  Pf  as  the  matrix  M  described  in  Lemma 


3.15 


Theorem 


3.16  shows  that  the  matrix  M  is  in  fact  a  five-diagonal  matrix;  and  Corollary  |3. 17  provides  the 


relation  between  the  eigenvectors  of  M  and  the  functions  i pn  defined  in 
Lemma  3.15.  Let  f  be  a  smooth  function  with  the  expansion  h  =  (ho,  hi,  ...)T  specified  in 


/(*)  =  J2hkPk(x)- 

k= 0 


(86) 


Suppose  that  ip  =  Dx(f),  with  the  expansion  r)  =  (Ho,  Hi,  ...)T  such  that 

OO 

vix)  =  '%2vkPk(x)- 


(87) 


k= 0 


Then, 


rj  =  Mh, 

where  the  matrix  elements  Mj ^  of  M  are  specified  via  the  formula 
Mjk  =  £pj(x)  ( Dx(Pt *))  (x)dx, 
with  0  <  j,  k  <  oo. 


Proof.  By  the  linearity  of  the  differential  operator  Dx  (defined  in  ( 69 ) ) , 

OO 

p(x)  =  ( Dx(f ))  (x)  =J2hk  iD^Pk))  (x)- 


k= 0 


Combining  (86)  and  (90), 


J2r,kP*(x)  =  J2hk  (Dx(P*))(x). 


k= 0 


k= 0 


Now,  by  multiplying  both  sides  of  (91)  by  P*  and  integrating,  we  have 


r»l  /  OO 


Vj  = 


5 Zhk(Dx(P *))  (x)  P/(x)dx, 


\k= 0 


(88) 


(89) 


(90) 


(91) 


(92) 
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and  by  linearity 


Vj 


oo 

£ 

k= 0 


hk  (£  P*(x)  ( Dx(p >*))  (®)d® 


(93) 

Q 


Theorem  3.16.  For  any  k  >  0, 


{DX(P*))  (x)  = 

— _ (k—l)2k2 _ “p* 

4v/2fc=3(2fc-l)v/2fc+T 


k3(l+/3) 


P>* 
7„ 


F2k-is/2k+l 

(— 4—6/3— 2fc/9(2+3£)+fc2(7+12(3+2/32)  +  (2fc3+fc4)(7+16/3+8/32))p5r, 

2(2fc— l)(2fc+3)  W 

(fc+l)3(l+/3) 


v/2FFTv/2fc+3  fc+! 

(fc+l)2(fc+2)2 
4v/2fc+T(2fc+3)%/2fc+5''  ^+2 


-ffc4-2(‘T)’ 


(94) 


in  (61). 

In  other  words,  M  is  the  five-diagonal  matrix 


where  Pj*  is  the  Normalized  Shifted  Legendre  Polynomial  defined  in  (T6j)  and  fi  =  is  defined 


Mk-2,k 
Mk-i  ,k 
Mk,k 
Mk+i,k 
Mk+2,k 


(k—l)2k2 

4s/2kNS(2k-l)^/2k+i 

k3(l+/3) 

v/2AETv/2fc+T 

(— 4—6/3— 2fc/3(2+3/3)+fc2(7+12/3+2/32)  +  (2fc3+fc4)(7+16/3+8/32)) 
2(2fc— l)(2fc+3) 

(fc+l)3(l+,3) 

V2k+IV2k+3 
(fc+l)2(fc+2)2 
4v/2FFT(2fc+3)v/2FF5  ‘ 


Proof.  By  the  definition  of  Dx  in  ( 69 ) , 


(95) 


(Dx(Pf))(x)  = 

=  ^  ((P  +  x)(fi  +  1  +  ®)®(1  -  x)-^P%(x)j  -2 x(x  + fi)Pf(x). 
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Using  the  chain  rule, 


(Dx(P*k))  (x)  = 


d 


=  ~r(P  +  x)(/3  +  l  +  x)  x(l  -x)-Pf(x) 


dx 


d 


d 


dx 


dx 

d 


+  (P  +  x)(P  +  l  +  x)—  x(l  -x)—Pk(x)  -2 x(x  +  P)P£(x)  = 


dx 


=(1  +  2x  +  2/3)  (x(l-x)-^P£(x)\ 

+  (x2  +  x(l  +  2/3)  +  (3  +  (32)  (x{l  -  x)-^Pf(x)^ 
-  2x(x  +  /3)Pfc(x). 


(97) 


Using  identities  (|13j),  (|14j)  and  ([15]) , 
(Dx(P*k))(x)  = 

(-1  +  k?k2PU(X) 


4(— 1  +  2/c)  (1  +  2fc) 

A^l  +  flP^s) 

1  +  2fc 

(— 4+7fc2+14fc3+7fc4— 6/3— 4fc/3+12fc2/3+32fc3/3+16fc4/9— 6fc/32+2fc2/32+16fc3/32+8fc4/32) 


2(  — 1  +  2/c)(3  +  2k) 


Pk(x) 


(1  +  fc)3(l  +  /3)P£+1(x) 
l  +  2k 

(l  +  k)2(2  +  k)2P*+2(x) 
4(1  +  27c)  (3  +  2k) 


(98) 

□ 


Finally,  substituting  ([6])  into  (98)  gives  (94). 

Corollary  3.17.  Suppose  that  hn  =  (/iq,/i", 
in  the  expansion  of  the  function  ipn(x)  defined  in  (63);  then,  hn  is  the  n  +  1-th  eigenvector  of 
M  : 


J 


is  the  vector  of  coefficients  defined  in  (If), 


Mhn  =  Xnhn, 


(99) 


where  M  is  the  five-diagonal  matrix  (95),  Xn  are  the  eigenvalues  of  the  differential  operator 
Dx,  and  k  =  0, 1,  2.... 


Proof.  By  (70),  ipn(x )  is  an  eigenfunction  of  Dx,  with  the  eigenvalue  Xn- 
{Dxflfn))  {x)  =  Xn^n{x), 


(100) 
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so  that 


dx  ( x  hkpkj  j  (*)  =  (dm)  (X) = xn  x  hkpk(x'y 


\k= 0 


Therefore,  by  Lemma  3.15,  we  obtain  (99) 


(101) 

□ 


3.9  A  relation  between  un,um,  and  the  ratio  an/ari 
Lemma  3.18.  For  any  n,  m  =  0, 1, ...., 

afn  =  Ig  u'n(t)um(t)dt 
an  /a  Un(t)u'm(t)dt 


(102) 


an  and  am  are  the  singular  values  (defined  in  (45))  of  Ca,b  (defined  in  (42)). 
Similarly, 


if  the  integrals  are  not  0;  where  un  and  um  are  right  singular  functions  (defined  in  (45))  and 
the  singular  values  (defii 

_  fo  /lp'n(x)lfm(x)dx 


an  fd  fin(x)fi'm(x)dx 


(103) 


if  the  integrals  are  not  0;  where  ifn  and  fim  are  defined  in  (63). 


Proof.  We  recall  from  (46)  that 


1  i  r°° 

un{t)  =  —  (C*(vn))  (t)  =  —  /  e~utvn(oj)du.  (104) 

JO 

Therefore,  the  derivative  of  un(t )  is 

i  r°° 

Unit)  =  —  /  i-u)e~utvniu])du.  (105) 

JO 

We  multiply  both  sides  of  the  expression  by  um(t),  integrate  both  sides,  and  change  the  order 
of  integration: 


’  1  / 

un(t)um(t)dt  =  —  / 
an  J  a 


i~co)e  wtun(w)do;  um(t)dt. 


By  rearranging  the  result,  we  obtain 

fb  otm 

/  u'n(t)um(t)dt  =  —  i~u)vniu;)vmiu)dw. 
J  a  an  Jo 

m ,  and  n  are  clearly  interchangeable,  so  that 


f 


Oi  i 

i~Uj)vniuj)vmiuj)duj  =  —  u'm(t)unit)dt. 


OLr 


(106) 


(107) 


(108) 


By  substituting  (|108)  into  (107),  we  obtain  (102).  Substituting  ([63])  into  (102)  we  obtain  (103) 


□ 
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4  Algorithms 


4.1  Computing  the  right  singular  function  un 

In  this  section  we  introduce  an  algorithm  for  the  numerical  evaluation  of  un(t),  the  n  +  1-th 
right  singular  function  (defined  in  (45))  of  Ca  &  (the  operator  defined  in  (42)). 


Step  1: 
Step  2: 
Step  3: 


Compute  hn,  the  n  +  1-th  eigenvector  of  the  matrix  At,  defined  in  (95) 


Compute  the  function  ipn(x)  from  hn,  using  the  expansion  specified  in  (75) 


Obtain  un(t)  from  ipn(x)  using  (63) 


Remark  4.1.  For  computations  in  precision  e,  the  vector  hn  =  (/iq  ,  /i”,  ...)T  is  truncated  at 
K,  such  that  \hj^\  <C  e  for  all  k  >  K.  By  lemma  3.13  the  coefficients  h £  decay  rapidly  as  k 


grows;  the  actual  position  of  the  last  significant  coefficient,  larger  in  magnitude  than  e,  is  given 
in  Figure  [6]  and  Table  [4]  in  Section  [5]  for  several  combinations  of  7  and  n. 

4.2  Computing  the  singular  value  an 

In  this  subsection  we  present  an  algorithm  for  computing  the  N  +  1  first  singular  values 


cuO)  07,  •••Oat  (defined  in  (45))  of  JZa,b  (the  operator  defined  in  (42)). 


Step  1:  Compute  the  first  singular  value  ao,  for  example,  by  the  formula 


oq  — 


fa  thu0(s)ds 

u0(t) 


(109) 


derived  from  (50),  where  t  G  [a,  b]  and  Uo(t)  is  computed  using  the  algorithm  from  the  previous 
section. 


Step  2:  For  every  n  >  0,  compute  an  from  an_i  using  the  relation  in  Lemma  3.18  and  the 
functions  un(t)  and  un-\{t)  computed  using  the  algorithm  in  the  previous  section. 
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5  Numerical  results 


In  this  section  we  present  results  of  several  numerical  experiments.  The  algorithms  for  comput¬ 
ing  the  right  singular  functions  un  and  singular  values  an  of  Ca ^  (the  operator  defined  in  (42)) 
were  implemented  in  FORTRAN  77,  using  double  precision  arithmetic,  and  compiled  using 
GFORTRAN. 

In  Figures  [lj  [2]  and  [3]  we  present  examples  of  right  singular  functions  of  the  operator  Ca,b, 
where  a  =  1  and  6  =  1.1,  b  =  10  and  b  =  100000  respectively. 

In  Figure  [4]  and  Table  [l]  we  present  the  singular  values  an  of  the  operator  £a,b,  for  several 
ratios  7  =  6/a;  an  depends  only  on  7  and  n  (see  remark  3.12).  In  table  [2]  we  present  several 
singular  values  smaller  than  IO-1000;  the  Fujitsu  compiler  with  quadruple  precision  was  used 
in  this  experiment. 

In  Figure  [5]  and  Table  [3]  we  present  the  eigenvalues  of  the  matrix  M  defined  in  (95). 

In  Figure  [6]  and  Table  [4]  we  present  for  several  combinations  of  7  and  n  the  position  of 
the  last  significant  coefficient  h £  in  the  expansion  defined  in  (75),  that  is  larger  in  magnitude 
than  e  =  10~16.  In  numerical  computations,  the  vectors  are  truncated  around  that  point  (see 
Remark  |4.1|) . 

In  figure  [5]  we  present  the  CPU  time  required  for  the  computation  of  the  expansion  of  the 
101-st  right  singular  function  u  100  of  £1,7,  for  varying  7;  The  experiment  was  performed  on  a 
ThinkPad  X230  laptop  with  Intel  Core  i7-3520  CPU  and  16GB  RAM. 
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Figure  1:  Right  Singular  functions  of  C\t m- 
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Figure  2:  Right  Singular  functions  of  £i,io- 


22 


Figure  3:  Right  Singular  functions  of  £1,100000  ■ 
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Figure  4:  Singular  values  an  of  Ca,b,  with  7  =  6/0. 
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Figure  5:  Magnitude  of  the  eigenvalues  of  the  matrix  M  defined  in  (95),  with  7  = 
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Last  significant  coefficient 


Figure  6:  The  position  of  the  last  significant  coefficient  larger  in  magnitude  than  10 


16 
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Figure  7:  CPU  time  required  for  computing  the  expansion  of  the  101-st  right  singular  function 
of  £1,7,  as  a  function  of  7.  The  experiment  was  performed  on  a  ThinkPad  X230  laptop  with 
an  Intel  Core  i7-3520  CPU  and  16GB  RAM. 
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Table  1:  Singular  values  an  of  Ca^ 
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Table  2:  Examples  of  singular  values  an  smaller  than  10 
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Table  3:  Eigenvalues  of  M  defined  in  |95 
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Table  4:  The  position  of  the  last  significant  coefficient  larger  in  magnitude  than  10 
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6  Conclusions  and  generalizations 


In  this  paper  we  have  introduced  effective  algorithms  for  the  evaluation  of  the  right  singular 
functions  and  singular  values  of  the  Truncated  Laplace  Transform  Ca,b- 

As  is  evident  from  Remark  2.3  and  the  more  detailed  discussion  in  m ,  the  right  singular 
functions  of  Ca are  an  efficient  basis  for  representing  decaying  exponentials  on  the  interval 
[a,  6], 

An  algorithm  for  the  computation  of  the  left  singular  functions  of  Ca  j,,  which  is  the  remain¬ 
ing  component  in  the  computation  of  the  SVD,  will  be  presented  in  a  future  paper.  Additional 
asymptotic  properties  of  the  Truncated  Laplace  Transform  and  of  the  associated  differential 
operators  will  also  be  discussed  in  a  future  paper. 
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